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We investigate two-species bosons in a two-dimensional square lattice by quantum Monte Carlo 
method. We show that the inter-species attraction and nearest-neighbor intra-species repulsion 
results in the pair-supersolid phase, where a diagonal solid order coexists with an off-diagonal 
pair-superfluid order. The quantum and thermal transitions out of the pair-supersohd phase are 
characterized. It is found that there is a direct first order transition from the pair-supersolid phase 
to the double-superfluid phase without an intermediate region. Furthermore, the melting of the pair- 
supersolid occurs in two steps. Upon heating, first the pair-superfluid is destroyed via a Kosterlitz- 
Thouless transition then the solid order melts via an Ising transition. 

PACS numbers: 67.80.kb 67.60.Bc 03.75.Mn 



Supersolid (SS) is an intriguing but counter-intuitive 
concept, for it is characterized by the co-existence of a 
diagonal long-range order (DLRO) of solid and an off- 
diagonal long-range order (ODLRO) of superfluid (SF). 
The possible existence of the SS phase has been stud- 
ied experimentally and theoretical since the early 1970's 
[l|-l3|. In recent years, numerical simulations have estab- 
lished the existence of the SS phase in a number of in- 
teracting lattice boson models [Jl. Recently, the concept 
of pair-supersolid (PSS) for two-species lattice bosons is 
proposed [5| . Here the pair-superfluid (PSF) replaces the 
SF to be the off-diagonal long-range order. In the PSF 
phase each species individually doesn't form the SF, in- 
stead the boson-boson composites form the PSF. It has 
been shown theoretically that PSF phase arises naturally 
in bosonic mixtures with inter-species attractions [6|-l8|. 
Some of the optimal candidates to realize the PSF phase 
include bosonic binary mixtures in optical lattices and 
dipolar particles in bilayer optical lattices. The later sys- 
tem is also an ideal candidate to realize the PSS phase 
because the dipole-dipole interaction can provide both an 
interlayer attraction and an intralayer nearest neighbor 
repulsion which give rise to the PSS phase [5|. 

In this work we consider the two-species Bose-Hubbard 
model in a two-dimensional (2D) square lattice, which is 
characterized by the Hamiltonian; 

^ = -^E ("I'^j + ^1 + ^■^) + ^ E <"J 

{ij) (uV 

+ yE<«-i) + ^E">'-^E<' w 

where a — a,b indicates the two species respectively, U 
{W) is the on-site intra-species (inter-species) interac- 
tion, and y > is the intra-species nearest neighbor 
repulsion. It is known that without the intra-species re- 
pulsion V, the Hamiltonian gives rise to PSF phase for 



T^ <0 
W > 




and super-counterfluid (SCF) phase for 



9g, in addition to the conventional Mott insu- 
lator and double-superfluid (2SF) phases. It has been 
proposed that the inclusion of intra-species repulsion V 
results in the checkerboard solid and PSS phases [5|. In 
Ref. 5] the existence of such a PSS phase is predicted by 
solving the time dependent Gutzwiller equations for the 
low-energy effective Hamiltonian of boson pairs. How- 
ever, the stability analysis of the PSS phase against phase 
separation is beyond such a treatment. Furthermore, the 
nature of the associated quantum phase transitions re- 
mains unknown and the precise phase boundaries are not 
determined. It is thus imperative to study the phase di- 
agram beyond Gutzwiller approximation and to investi- 
gate the nature of associated quantum phase transitions. 
Additionally, so far the theoretical investigations of the 
PSS phase have focused on it's existence and the stability 
at zero temperature. Since it is known that typically the 
thermal melting of a SS state contains two distinct tran- 
sitions ll|, |l2| , a precise characterization of the thermal 
transition out of the PSS state is also called for. 

In this work we perform large-scale quantum Monte- 
Carlo (QMC) simulations to study the quantum and 
thermal phase transitions of the two-species Bose- 
Hubbard model characterized by Eq. ([T|). In particu- 
lar, we employ a multi-worm algorithm which is similar 
to the method proposed in Ref. [13|. We shall briefly 
summarize our main results before we discuss them in 
detail: We obtain the ground state phase diagram in the 
/i, t plane as shown in Fig. [T] Two checkerboard solid 
lobes with density n^ = 1/2 and n^ = 1 respectively are 
identified. Here ria = X]j("f)/^^ i^ the average parti- 
cle density of a-species boson and L is the system size. 
We find that the PSS phase appears between these two 
lobes while the PSF phase appears between the n^ — 1/2 
lobe and vacuum. We observe that all transitions into 
the 2SF phase are first-order in the parameter region we 
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FIG. 1. (Color online) Ground state phase diagram for 
the two-species Bose-Hubbard model for W — —0.95(7 and 
V = 0.025(7. First-order (continuous) transitions are shown 
by filled (open) symbols and are estimated using data with 
L — 30 {L — 24). The black dotted lines indicate the phase 
boundaries determined by the mean-field solution as described 
in Ref . [a] . The density profile of the checkerboard solid lobes 
are sketched inside the lobe. 



FIG. 2. (Color online) Finite temperature phase diagram 
for the two-species Bose-Hubbard model for W — —0.95U, 
V = 0.025(7, and /i = —0.435(7. First-order (Ising) transition 
is shown by filled triangle (open square), while KT transi- 
tions are shown by open circles. Finite-size scaling is used to 
determine the critical temperature as described in the text. 
Black dotted lines represent the estimated phase boundaries 
near the multicritical point. Lines are guides to the eyes. 



Studied. These include the PSS-2SF, PSF-2SF, solid- 
2SF, and vacuuni-2SF transitions. In contrast, the PSS- 
solid and PSF-vacuum transitions are continuous. But 
the PSF to Ha = 1/2 solid transition is first-order. Fur- 
thermore, we establish that the transition from PSS to 
2SF is direct without an intermediate region. We also 
obtain the finite temperature phase diagram on the t, T 
plane for a fixed fj, — —0.435 as shown in Fig. [21 We find 
that the melting of PSS phase contains two distinct tran- 
sitions: Upon heating first the PSF is destroyed via a KT 
transition. When heated further the solid order melts via 
an Ising transition. Finally, the finite-temperature solid 
to 2SF transition is first order. 

In the following we discuss our simulations and findings 
in more detail. In all calculations we set (7 = 1 as our 
energy unit. We note that the system can exhibit three 
kinds of ODLRO corresponding to three kinds of super- 
fiuid density. It is known that single species SF density 
pj can be estimated by p^ — ((Wcr) )/(4/3i), where W^r 
is the winding number for a species boson and /3 is the 
inverse temperature [14]. The PSF (SCF) density can be 
identified via the sum (difference) of the winding num- 
ber with the formula Ps — ((W±) )/(2/3i), where 
W± = (Wa ± Wb)/2. The PSF phase is identified by a 
finite PSF density {p^^^ ^ 0) together with a zero SCF 
density (pf^^ = 0), which characterizes a perfectly corre- 
lation between bosons of different species. The system is 
in a 2SF phase if all three kinds of superfluid densities are 
nonzero. For small enough t, all kinds of superfluid den- 
sities are zero and the ground state may acquire a DLRO, 



which can be detected by measuring the structure factor 



Due to the symmetry of the model one has Sa ~ Sb and 
the index a can be dropped for simplicity. The system 
is in the PSS phase if the structure factor is nonzero 
while the system also has perfectly correlated superfluid- 
ity (pPSF ^ and pfCF ^ 0). 

When U + W,V <^ U, a low-energy effective Hamilto- 
nian of pairs can be derived using second order pertur- 
bation theory. Within the mean-field theory the paired 
order parameter can be determined by solving the self- 
consistency equation [5] . When U + W < AV the mean- 
filed analysis predicts that both the Ua = 1/2 and Uc = 1 
insulating lobes have checkerboard order and PSS/PSF 
phases may exist outside the lobes. To firmly estab- 
lish the existence of the PSS phase and to compare pre- 
cisely with the results in Ref. [5] we set W — —0.95 and 
V — 0.025 throughout this work. As shown in Fig. [1] we 
find that mean-field results substantially deviate from our 
QMC results and mean-filed solutions overestimate the 
region of solid phase by a large margin. The PSS/PSF 
regions identified by QMC are smaller than the results 
by time-dependent Gutzwiller equations and no PSS re- 
gion is found below the Ua = 1/2 lobe. In addition for 
the first lobe we observe a very weak reentrant behavior, 
which is similar to what has been found in ID, but such 
a reentrant behavior is absent for higher lobes [l0[ . 

Now focusing on the quantum transitions out of the 
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FIG. 3. (Color online) (a) Energy per site E/N, (b) structure 
factor S/N, (c) PSF density p^^, and (d) SCF density p?^^ 
for p = -0A35U, /3 = 1000, and L = 24. The hysteresis ef- 
fects were obtained by increasing (decreasing) t in small steps 
starting from a configuration inside the PSS (2SF) phase. The 
last configuration generated in the run at a particular t are 
used as the starting configuration for the next value of t. The 
critical hopping tc is extracted from the crossing of the energy 
branches. The QMC error bars are smaller than the symbols. 



FIG. 4. (Color online) PSF density in the vicinity of the 
KT transition temperature T^t for t = 0.045 with system 
sizes 8 < L < 32 and the infinite system (solid line). Inset: 
K(r) in Eq. ((3]) for different pairs of system sizes (Li,Z/2) = 
(8, 16)(blue circle), (8,24)(green squre), (8,32)(red triangle) 
and (f 6, 32)(magenta diamond). The dash line shows the liner 
fitting for k(T) = l + c{T^^^ -T) with T^f^ ~ 0.00283(7) and 
c ~ 134.88(1) respectively. The QMC error bars are smaller 
than the symbols. 



PSS phase. Similar to the continuous PSF-vacuum tran- 
sition, we find that the PSS-solid transition is also contin- 
uous. The PSS-2SF transition, however, is strongly first 
order. Near the first order phase boundary the Marko- 
vian dynamics of the simulation suffers from the critical 
slowing down due to the large tunneling time between 
the competing phases. Interestingly, it is pointed out 
recently that such a large tunneling time can be used to 
locate the phase boundary accurately by determining the 
average energy of the two different phases in the proper 
manner [15| . We first determine the energy of the PSS 
phase by starting the simulation at a value of t inside 
the PSS phase and then increase t in small steps to go 
to the 2SF phase. During this procedure we always use 
the last configuration generated in the run at a particu- 
lar t as the starting configuration for the next value of 
t. For large enough system the tunneling time is longer 
than the simulation time and the procedure ensures that 
the system stays in the metastable PSS state even when 
the phase boundary has been crossed. By starting from a 
value of t inside the 2SF phase and decrease t in the same 
fashion, the energy of the 2SF phase is obtained. The 
crossing point of the two energy branches then gives tc- 
In Fig. [Slja) we show the energy per site for /i — —0.435, 
/3 — 1000, and L — 24. The crossing point of two en- 
ergy branches gives an estimation of tc — 0.0477(2). In 
Fig.[2Ib), (c), (d) we also show the measurements of S/N, 
p^^^ , and p^'-^^ respectively, where the hysteresis behav- 
ior is clearly observed. 

We now turn to the finite temperature properties. For 



PSS phase, it is expected that the solid order can per- 
sist up to some finite temperature T^ beyond which the 
solid order melts due to thermal fluctuation. On the 
other hand, while in 2D one cannot break the continuous 
symmetry at finite temperature, quasi-long-range order 
can persist up to a Kosterlitz-Thouless (KT) transition 
temperatures. Consequently the PSF density p^^^ re- 
mains finite up to the critical temperature ^xf^ and a 
persistence of the PSS state at finite temperature is ex- 
pected. However, it is not straightforward to see whether 
the solid order and the PSF density disappear simultane- 
ously, or which order will be destroyed earlier by thermal 
fiuctuation. It is also not clear whether an intermedi- 
ate SS or PSF region can appear at finite temperature. 
Furthermore, it is expected that for 2SF phase the su- 
perfiuid density p^ should also remain finite up to the 
critical temperature T^j,, but it is not clear whether T^rp 
will connect smoothly to T^^F as one varies the system 
parameters across the PSS-2SF phase boundary. 

To accurately determine the KT transition tempera- 
ture we utilize the KT renormalization group equations. 
It is known that for an infinite system there should be a 
universal jump at the transition temperature. For finite- 
size system, however, the transition is smeared by loga- 
rithmic finite-size effects. Here we follow the proposed 
method in Ref. [ll[ to determine T^^ (Trt^) in the 
thermodynamic limit. We define R1{L,T) = tTTp1{L)/T 
{R^^'^{L,T) = t7rpPSF(^)/2r) and study the finite-size 
scaling using KT renormalization group equations in the 
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FIG. 5. (Color online) Rescaled structure factor 

S{ii,Ti)L'^^'^ /N versus temperature for various system sizes 
near the solid melting temperature. The crossing point indi- 
cates the critical temperature Tc — 0.0335(5). Inset: Data 
collapse with rescaled temperature (T — Tc)L^^'' . The Ising 
critical exponents /3 — 1/8 and u = 1 are used. The QMC 
error bars are smaller than the symbols. 



destroyed before the melting of the solid order occurs. We 
note that PSF and SF density disappear simultaneously 
while the solid order is still very strong above T^f^ . This 
indicates that there is a direct PSS-Solid transition with- 
out an intermediate SS phase. We also observe that the 
finite-temperature solid-2SF transition is a direct first or- 
der transition without an intermediate region. The first 
order nature is expected due to the different symmetries 
in different phases. 

Finally we study the melting of the sold phase. It 
is expected that at the critical temperature Tc the solid 
melts via 2D Ising transition. To confirm the Ising nature 
of the transition we use the finite size scaling hypothesis 
for the structure factor S/N = L-^I^/''T{{T - Tc)L^'''). 
The Ising universality class can be checked via the data 
collapse for {S/N)L^I^/'' v.s. (T - Tc)L^''' for various L 
with the Ising exponent /3 = 1/8 and v = 1. In Fig. [5] 
we plot the re-scaled structure factor as as a function of 
the temperature, where the crossing point provides an 
estimation of the critical temperature Tc- Data collapse 
is clearly observed if the temperature is also rescaled, as 
shown in the inset of Fig. [5] The Ising nature of the 
transition is hence clearly confirmed. 



integral form [if 



41n(L2/ii 



R{Li,T) 



dt 



B.(L,X) t^(ln{t)-K{T))+f 



(3) 



Here k{T) is a system size independent parameter, which 
is analytic in terms of temperature. It is expected that 
k{T) ~ l-t-c(rKT-T) for T < Tkt- For each pair of sys- 
tem sizes, one can determine a curve for k{T). The data 
supports a transition of the KT type if k(T) obtained 
from different pairs of the system sizes collapse into a 
straight line around n — 1. The transition temperature 
is then determined by k(Tkt) — 1 and the superfluid 
density in thermodynamic limit can be determined by 
the equation 1/R + Ini? = k{T). 

In Fig. |4] we show the PSF density as a function of 
temperature for different system sizes for t = 0.045 and 
fi — —0.435. In the inset we show k(T) obtained from 
different pairs of system sizes. The data collapse and the 
smooth analytic behavior of k{T) confirm the KT na- 
ture of the transition, and the transition temperature is 
found to be T||P ~ 0.00283(7). The PSF density in ther- 
modynamic limit with universal jump is plotted as solid 
line. Similar analysis were applied to the 2SF-normal 
fiuid transition. By identifying T^rp and Tj^f^ at differ- 
ent t we find that in general T^rp is much higher than the 
T^jF J as shown in Fig. [2] They don't smoothly connect 
to each other as t is varied across the zero temperature 
PSS-2SF phase boundary 

In Fig. [2] we plot the finite temperature phase diagram 
on t—T plane for the cut at fi = —0.435. We find that the 
melting of PSS occurs in two distinct steps: the PSF is 



SUMMARY 

In summary we have studied the two-species Bose- 
Hubbard model with on-site intra-species (inter-species) 
repulsion (attraction) and nearest neighbor repulsion in 
a 2D square lattice. We obtain the ground state and 
finite temperature phase diagrams and firmly establish 
the existence of PSF and PSS phases within this model. 
We accurately determine the phase boundaries and the 
nature of the phase transitions. Interestingly, we find 
that the transition from PSS to 2SF or solid phase is 
direct without an intermediate phase. Furthermore, the 
melting of the PSS occurs in two steps: the PSF is first 
destroyed via a KT transition then the solid order melts 
via an Ising transition. 
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